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ABSTRACT 

Using high resolution adaptive mesh refinement simulations in 3D, we investigate the 
formation of relativistic jets from rotating magnetospheres. Here, we focus on the 
development of non-axisymmetric modes due to internal and external perturbations 
to the jet. These originate either from injection of perturbations with the flow or from 
a clumpy external medium. In the helical field geometry of the accelerating jet, the 
m=l to m=5 modes are analyzed and found to saturate at a height of ~ 20 inner 
disk radii. We also discuss a means to control artificial amplification of m = 4 noise 
in the cartesian simulation geometry. Strong perturbations due to an in-homogeneous 
ambient medium lead to flow configurations with increased magnetic pitch and thus 
indicate a self- stabilization of the jet formation mechanism. 

Key words: accretion, accretion disks - ISM: jets and outflows - MHD - galaxies: 
active - galaxies: jets - relativity 



1 INTRODUCTION 

Relativistic jets from active galactic nuclei (AGN) are host 
to a multitude of physical processes that lead to emission 
of radiation with the highest energies in the universe (e.g. 
Bottcher 2007 ) . The observed non-thermal emission of high- 
energy particles is subject to beaming in a highly collimated 
jet of plasma moving with Lorentz factors in the range 
10-40 ( jJorstad et al.||2005| |Lister et aLp009) . The lead- 
ing paradigm of jet formation suggests that the accelera- 
tion of the bulk flow can be understood in the relativistic 
magnetohydrodynamic (RMHD) approximation, where the 
jet is endowed with a large-scale current circuit. Current 
axisymmetric simulation models of jet forma tion (|Komis- 
sarov et al.] [20071 |Tchekhovskoy et al.||2008| |Porth et af 



201 1| have matured to an overall agreement on the underly- 



ing process of jet formation. However, since nature does not 
obey axisymmetry, 2D modeling falls short in two important 
aspects: 1. The axisymmetric induction equation (3$ = 0) 
lacks a mechanism to transform toroidal magnetic field back 
into poloidal field. Therefore the toroidal field can only be 
amplified, leading to a potential overproduction of B<j> and 
the collimating pinch force. This might cast a shadow of 
doubt on the validity of results from axisymmetric jet for- 
mation simulations. 2. In order to also study the stability 
of jets, a three dimensional treatment is necessary. It is well 
known that instability of current carrying plasmas is con- 
nected to non-axisymmetric perturbations of the toroidal 
magnetic field ( |Bateman|1978 ). Among these current driven 



instabilities, the m — 1 mode known as the "kink" is the 
most violent one. It leads to a helical displacement of the 
flow from the axis of the plasma-cylinder. In the context of 
young stellar objects, the helical kink can yield an explana- 
tion for the wiggly structure observed in some stellar jets 
(e.g. HH 46/47) ( |Todo et al.(l993| |Lery eTal^OO^ . If not 
brought to a halt by regulating non- linear mechanisms, the 
exponential instability growth must lead to complete dis- 
ruption of the jet. At the presence of instability, jets can 



dissipate magnetic energy via reconnection (e.g. Drenkhahn 
\k Spruit|2002||Lyutikov fe Uzdensky|2003t and shocks (e.g. 



Blandford &; Konigl ( 1979 )), leading to heating, acceleration 
of high-energy particles and radiation. As noted by |Heinz fc| 
Begelmanl (120001) andlDrenkhahn & Spruit I (120021), the dissi- 



pation into a fully tangled magnetic field could also promote 
efficient quasi-thermal acceleration of the bulk flow out of 
magnetic enthalpy. In principle, the current-driven instabil- 
ities can be accompanied by other types of instabilities such 
as the Kelvin-Helmholtz instability (KH) caused by shear 
between jet and ambient material. The KH instability can 
thus lead to an efficient mixing of the jet and environment. 
However, in the presence of strong magnetic fields, growth of 
the KH modes is strongly suppressed ( |Keppens et al.|1 999 ) . 
In particular, toroidal fields appear to hinder mixing and 
thus exert a stabilizing influence ( |Appl &; Camenzind| 1992 
Mignone et~aL|2010 ) 



As pointed out by McKinney & Blandford ( 2009 ) , when 



the well-known Kruskal-Shafranov (KS) instability criterion 
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for cylindrical force-free equilibria is applied to relativistic 
jets, we obtain the result that jets become unstable already 
at the Alfven point za — 10rs (where > B p and z > r; 
rs denotes the Schwarzschild radius) - before accelerating to 
highly relativistic velocities. This is in stark contrast to the 
finding of some AGN FR-II jets propagating unperturbed 
out to distances of 10 7 rs. 

Linear stability analysis in non-relativistic, but oth- 
erwise complete MHD was conducted for current-free and 
current-carrying jets with helical magnetic fields by |Appl"fc] 
Camenzind| p992 ) . The stability of relativistic MHD jets on 
the other hand is still not fully addressed in linear analysis 
and subject of ongoing research. Considering linear analy- 
sis of force- free cylindrical configurations, Ist omin &: Pariev] 
(1994) and Istomin & Pariev (11996) concluded stability of 



B z — const, jets with respect to axisymmetric and helical 
perturbations. Begelman (1998) on the other hand demon- 
strated the violent instability of the m — 1 mode which 
he postulated as a possible solution to the long-standing 
a— problem. Also |Lyubarskii| ( |1999[ ) stressed the importance 
of the m — 1 mode by considering more realistic configura- 
tions with decreasing flux. However, Lyubarskii] (fl999 ) and 
recently Narayan et al. (2009) note the time-dilated slow 
growth rate of the kink. This can lead to a substantial dis- 
placement between launching and dissipation scales and thus 
yield a possible explanation for the extent of the blazar zone 

of AGN. 

Tomimatsu et al. ( |2001 ) could extend the classical KS 



criterion and demonstrated a stabilizing effect of the rela- 
tivistic field line rotation. Their analysis (TKS) yields the 
simple criterion for instability 
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where Q denotes the angular velocity of the field line Q — 
(v<j> — VpBff,/ B p )/r, such that the toroidal field strength also 
has to overcome the stabilizing electric field. The asymp- 
totic relation for relativistic jets B^ ~ —rQ/cB p suggests 
marginal stability of the unperturbed flowQ Taken at face 
value, the TKS criterion thus suggests that relativistic jets 
are always on the verge of instability with the ultimate fate 
strongly depending on details of the modeling. In more re- 
cent studies, several authors investigated the stabilizing in- 
fluence of environmental effects such as shear (e.g. |Mizuno| 



et al. 2007), external wind with relativistic bulk motion 



^Hardee & Hughes 2003) and sideways expansion (Rosen & 



Hardee 2000 ) , emphasizing the influence of modeling details 



for jet and ambient material. While there is now growing 
consensus that the kink instability can also operate in the 
relativistic regime, to answer whether it can grow indefi- 
nitely to finally disrupt the jet or rather saturate, requires 
a study of the non-linear evolution via numerical simula- 



tions (e.g. Mizuno et al. 2009 Mignone et al. 2010 Mizuno 



et al.||2011| ). Also |Q'NeiU et al.|fl2012D studied the stability 
of various local jet configurations in force-free, pressure con- 
fined and rotational equilibriua. These authors found pres- 
sure confined models corresponding to astrophysical jets fur- 
thest from the origin to be the most unstable ones. 

1 However, this result should be handled with care, since the 
TKS criterion strictly only applies in the sub-Alfvenic region of 
the flow; see also the discussion by McKinney & Blandford ( 2009| ). 



On the other hand we would like to point out that jet 
stability can best be studied by including the initial ac- 
celeration and collimation region of the jet. Thereby, self- 
consistent helical magnetic fields are obtained and the num- 
ber of ad-hoc assumptions for the jet base can be reduced 
to a minimum of physically well motivated choices. |McK-| 
|inney &; Blandford| ( |2009| ) were the first to present global 
3D simulations of jet formation including a turbulent accre- 
tion disk. In their seminal paper, no significant disruption 
or dissipation out to scales of 10 3 rs was observed, however 
by using comparatively low resolution GRMHD simulations 
with 256 x 128 x 32 cells in a spherical r x x grid. Al- 
though launching jets directly from the turbulent accretion 
flow promises high realism, a systematic study of jet insta- 
bility growth can hardly be performed this way. For non- 
relativistic disk jet s, |Ouyed et ah] (|2003|, |Anderson et al 



(2006), Moll et al. (2008) and Staff et al. (2010) followed an 



alternative approach by treating the rotating magnetosphere 
as a fixed-in time injection boundary. As a natural general- 
ization of previous axisymmetric work, the boundary can be 
used to model the corona of a Keplerian accretion disk. Also 
we will follow this strategy for the case of relativistic jets, 
as it allows to systematically control both fixed in time and 
time variable injection conditions. 

In this paper, we show our first results concerning the 
stability of relativistic jets near the launching region. We will 
discuss 3D RMHD simulations exposed to non-axisymmetric 
perturbations triggered by the accretion disk and due to jet- 
cloud interactions. To complement previous considerations 
of axisymmetric (2.5D) jet formation outlined in |Porth fe| 



Fendt| ([20Toj) (PF I) and |Porth et ah] ( poTT) (PF II), we 



focus mainly on non-axisymmetric features within the jet. 



2 MODEL SETUP 

The jet is launched from a rotating inlet resembling the 
corona of an accretion disk similar to PF I. The initially 
purely poloidal magnetic field is transformed into a helical 
shape giving rise to a global electric current system which 
accelerates and collimates the flow into a jet. Since the disk 
enters our description essentially as a boundary (allowing for 
outgoing Alfven and fast waves) we do not investigate the 
back reaction from the jet to the disk. Due to the increased 
complexity of the three dimensional treatment compared to 
previous 2.5D work, we neglect stratification caused by a 
gravitational source term, since the flow is accelerated above 
escape velocity very near the disk and sonic surface (see also 
PF I). Note that in the following, the axis of symmetry is 
denoted as the cartesian 'y- ax i s ' (not the 'z-axis' as usual). 



2.1 Initial and Boundary Conditions 

We initialize the domain with constant values for density 
and pressure p = po = 0.01, p = po, threaded by an ini- 
tial poloidal magnetic field of monopolar shape with the 



2 Although a fourth order limiter was used, we note here that 
an accurate study of the higher-order mode growth in the jet 
requires sufficient actual angular resolution (McKinney Sz Bland-| 
|ford| ( [2009} report 20% convergence in the m = 1, 2, 3 power when 
comparing 16 and 32 angular cells). 
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Figure 1. Sketch of the computational domain in the x-y plane 
(left) and the x-z plane (right). Within r ou t, rotating injection 
conditions are applied. The monopolar magnetic field is indicated 
by dashed lines in the xy plane and has its "source" outside of 
the domain at (x, y, z) = (0, —4, 0). 



"source" outside of the domain at y = — 4. The initial pres- 
sure follows from the plasma /3 parameter at the inner disk 
edge r = 1 to po = f3/2Bi, where B\ denotes the corre- 
sponding poloidal magnetic field strength at (r,y) = (1,0) 
and we set to B\ — 1 for convenience. A sketch of the com- 
putational domain is shown in figure [I] 

To avoid the poorly determined states in the corners 
of the domain where disk-boundary and outflow boundary 
would meet, we truncate the disk inlet at r ou t- Not trun- 
cating the disk introduces additional m = 4 noise triggered 
predominantly at the corners of the domain. At the injec- 
tion boundary within r < r ou t we assign boundary condi- 
tions for a rotating magnetosphere, injecting a sonic flow 
aligned with the poloidal field- lines. To provide the axisym- 
metric boundary constraints at the bottom plane, we first 
transform to cylindrical coordinates around the jet axis in 
the middle of the cartesian domain, compute the required 
quantities and subsequently transform back to the carte- 
sian domain to assign the updated boundary conditions. We 
specify fixed-in time axisymmetric profiles for the five quan- 
tities p d = Po,Pd = lOO/Oo = l,V(i, = ru(r),v p = c s ,^ = 
to account for the downstream Alfven and fast magnetosonic 
waves leaving the domain across the boundary. In addition, 
the magnetic flux B y is fixed to the initial profile. Since 
the evolution of B y is already suppressed via the choice of 
E<f, = 0, this does not represent an additional boundary con- 
straint. The radial field component is obtained by a field ex- 
trapolation satisfying = of the domain values B r , given 
the fixed profile of the vertical field component B y . The 
toroidal field B<j> follows from the domain via zero gradient 
extrapolation to yield also j r = 0. This way, the boundary 
current is under control and no spurious current sheets are 
allowed to arise at the bottom plane injection boundary. For 
the rotation profile, we adopt 



uj(t) = 0.5c/Vh 



1 



-3/2 



r < r in 

1 < r < rout 



(3) 



representing an inner solid-body rotation law and an outer 
Keplerian profile within n n = 1 as in Tchekhovs koy et al.| 
( 2008 ) . Beyond r ou t , we simply freeze the initial conditions 



in the boundary. 

We should note here that in a general relativistic treat- 



ment, frame dragging of a monopolar magnetosphere leads 
to an approximate solid-body field-line rotation with Q c± 
0.5 Qh, where Qh denotes the outer horizon angular veloc- 
ity ( |Blandford fe Znajek|[l977| |Komissarov"||2001| ). In the 
extreme Kerr case, this would yield Q ~ 0.25 c/n n hence we 
are even exaggerating the central field line rotation and thus 
the putative Blandford & Znajek power output. Naively, one 
could thus expect a fast axial jet to arise in our simula- 
tions. However, since the mass influx into the domain is in 
fact increasing as one approaches the origin (given constant 
poloidal disk- velocity and density) and due to the neces- 
sity of vanishing Poynting flux on the axis, our simulations 
will not produce a high Lorentz factor spine. This would re- 
quire an accurate modeling of the mass flux in the plunging 
region using a GRMHD approach which we must defer to 
future work. 

For the set of parameters shown here, /3 — 0.01, V0, r =i = 
0.5c and pd = 1, the injection speed follows as c s ~ 0.1c and 
hence the maximal |Michel (1969) magnetization parameter 
can be written 



0~M 



47rTpdC s c 



3 * 



(4) 



It becomes ctm = 2.5 giving an upper limit of the Lorentz 
factor to be obtained by the flow. Higher values of ctm can 
be achieved by reducing the disk density or the plasma- f3 as 
the magnetization depends on these parameters according 
to cr M oc p- 1/2 /3~ 1/2 . 

Due to the accumulation of numerical errors and de- 
creased resolution compared to previous axisymmteric stud- 
ies, steep gradients of the rotation law at the inner disk edge 
must be avoided if no special treatment for the velocity pro- 
file is given|^]lt is also for similar stability reasons that we 
specify w as a boundary constraint instead of the field rota- 
tion law Q as customary. 

To close the set of ideal special relativistic magnetohy- 
drodynamic equations (RMHD), we employ the equation of 
state proposed by Meliani et al.|(|2004|, approximating an 



ideal one-component monoatomic Synge (1957) gas (see also 
equation 6 of |Keppens et al. (2012)). 

Additional passive tracer scalars i\ — L3 are advected 
with the flow for the purposes of refinement and post- 
processing. 



2.2 Numerical Grid Setup 

We perform simulations in a cartesian adaptive grid using 



MPI-AMRVAC QKeppens et al.||2012| . With the adaptive 
mesh refinement in a cartesian domain, a uniformly high 
resolution for all regions of interest is obtained without the 
directional biases present in stretched grid simulations (e.g. 
Moll et aL]|2008| [Mignone et al.||2010|, especially, a high 



resolution can be achieved also at the jet head which will 
prove necessary to resolve the helically displaced tip of the 
jet. The cartesian discretization and quadrant al symmetry 
of the grid can however introduce notable noise leading to 
a "pumping" of the m — 4 mode. We comment on the mea- 
sures taken to analyze and circumvent such spurious effects 
further in section [3^21 



See also the discussion in Appendix 3 of |Ouyed etaT] ( |2003| ) 
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The largest domain size considered in the following ex- 
tends over x e [-32,32], y e [0,128] and z e [-32,32] 
in units of the inner disk radius. A base resolution of 
n x xn y xn z = 48 x 96 x 48 cells is chosen, adapt ively refined 
by four additional grid levels (three for the smaller domain 
case). Thus 12 grid cells per inner disk radius are achieved, 
totaling in 240 cells across the entire jet inlet radius. To our 
knowledge this high resolution is unprecedented in 3D jet 
formation simulations. The effective resolution for the large 
domain is 768 x 1536 x 768 = 9.06 x 10 8 cells and we observe 
a grid filling with ~ 3 x 10 8 cells at the termination time. 



2.2.1 Refinement Strategy 

Refinement to the highest level is enforced for the disk inlet 
r < Tout;?/ < 1 and for the region around the axis x, z < 
1; y < 10 to resolve the steep gradients of the monopole field. 
The jet is refined based on the scheme proposed by [Lohner 
( 1987 ) with weighted contributions of density, Lorentz factor 



and magnetic field strength. 

As additional criterion, the initial domain is coarsened 
to the lowest level until the jet comes within reach. This de- 
cision is based on the Lorentz factor (T > 1.05) and the pas- 
sive tracer scalars described earlier. The criterion in terms 
of the Lorentz-factor is required in order to avoid acciden- 
tal coarsening of the torsional Alfven wave launched at the 
switch on of disk-rotation; it does not transport jet-inlet ma- 
terial and is thus unaffected by the tracers injected along the 
jet. With this strategy, the fast jet is always resolved appro- 
priately, while the enforced coarsening of the initial domain 
significantly speeds up the simulations. Especially when the 
domain is initialized with a non-homogeneous medium, the 
coarsening-strategy is crucial. 



2.3 Perturbations 

In order to investigate the behavior upon instability and 
to break the quandrantal symmetry of the grid, the ini- 
tially axisymmetric setup needs to be distorted by non- 
axisymmetric perturbations. The physical origin of the per- 
turbation can be found in a non-axisymmetric evolution of 
the accretion disk, for example due to orbit of vortices and 
quasi-periodic oscillations 



van der Klis et al. 1985), or in 



a non-homogeneous external medium. In active galactic nu- 
clei, the presence of such a medium within the central region 
is well established as the broad line region that is possibly 
comprised of clouds orbiting the black hole at high veloc- 
ity (~ 0.1c) (e.j 



Davidson & Netzer 1979 Araudo et al 



2010 ) or in manifestation of a clumpy torus surrounding the 



accretion disk (e.g. Dullemond & van Bemmel 2005| . The 
eventual collision of such clouds with a relativistic jet was 
already proposed by Blandford &; Konigl (1979) as a mecha- 
nism to explain transient features in compact radio sources. 
The scenario of jet-cloud collision is thus of interest not only 
in respect to jet stability, but also concerning the triggering 
of particle acceleration at the shock surface. We take the 
presence of a clumpy ambient medium as a motivation for 
our second setup of jet perturbations, where we model how 
the emerging jet propagates through such a density struc- 
ture. 



Table 1. Parameter summary of the 3D simulations. Simulation 
names indicate box size (L,M) and the nature of perturbation (m 
- mode injection, d - density perturbations). 



L3D 

L3Dm 

M3D 

M3Dd 

M3Dmd 



0.01 
0.01 
0.01 
0.01 
0.01 





0.02 




64 X 128 X 64 
64 X 128 X 64 
32 X 64 X 32 
32 X 64 X 32 
32 X 64 X 32 



168 
168 
103 
189 
222 



2.3.1 Mode Injection 

To model non-axisymmetric features of the accretion disk, 
we employ perturbations to the rotation velocity uj + Auj 
following a mode decomposition. This method was already 



successfully applied by Rossi et al. (2008| and |Mignone et al 



2010J) . In particular, we adopt 



Au(r,t) 



32 



3 8 

cos(ra0 + u(l)uo(r)t + (f>o(l)) 

m=0 1=1 



(5) 



to obtain modes with sub- and super-Keplerian frequencies 
u(n) e {0.5, 1, 2, 3, 0.03, 0.06, 0.12, 0.25} featuring a random 
phase offset <po(l). The maximum amplitude of the pertur- 
bation is set to €uj = 2%. 



2.3.2 Clumpy medium 

A static clumpy medium is modeled as density perturba- 
tion by prescribing the size spectrum in Fourier space and 
subsequent transformation of the random-phased Fourier co- 
efficients to real space. To obtain a particular cloud size, the 
following spectrum is adopted: 



/(*0 = K 



_i g_+i 

2s 



k s 

n/ mi 



k ^ &min 



(6) 



where k m in relates to the cloud size A and domain size via 
A = max(Ax, Ay, Az)/k m in- The highest frequency /c max is 
limited by the base-level resolution of the grid to k max < 
0.5min(n cc , n y , n z ) as we do not allow refinement on the ini- 
tial condition. We adopt a slope of the cloud size function 
s = 5/3. In the following, the maximal cloud density is set 
to 100 times the ambient density po. 



3 RESULTS AND DISCUSSION 

We now describe the simulations performed and the analysis 
of non-axisymmetric features. The non-linear temporal and 
spatial evolution of angular modes is discussed and we show 
evidence for jet self-stabilization in the launching region. 



3.1 Overview of the Simulations 

Table [l] summarizes our parameter runs and Figure [2] shows 
a rendering of our reference unperturbed solution labeled as 
M3D. 

Various flow quantities in the z — plane containing 
the jet "axis" are shown in Figure [3] Here, we also illus- 
trate the momentary run of the critical surfaces as in PF 
I. Note that these are strictly only valid under the assump- 
tion of stationarity and axisymmetry. Compared to PF I, 
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Figure 2. Rendering of the reference simulation, showing cut 
open pressure isocontour colored according to Lorentz factor with 
magnetic field lines (run M3D). 



our new simulations are not evolved long enough to adopt 
the near- stationary state and still reflect the initial transient 
phase of the sudden spin up of the disk. None the less, the 
flow transcends the critical Alfven surface (blue contour) 
and becomes super-luminal (black contour) but stays sub 
magnetofast (red contour) within the domain. 

At this scale, the Lorentz factor of the disk wind is 
still moderate (r < 2), it rises above T ~ 2 only in regions 
localized adjacent to the high pressure backbone. 



3.2 Mode Analysis 

We now first evaluate the impact of artificial m = 4 pump- 
ing due to the quadrantal symmetry and quantify the growth 
of non-axisymmetric modes within the jet formation region. 
To first give an impression of the azimuthal variations, a 
qualitative comparison of the unperturbed simulation L3D 
with run L3Dm is shown in slices of selected quantities mid- 
way across the jet in Figures [4] and [5] When no measure of 
perturbing the quadrantal symmetry is taken, we observe 
multiples of the m — 4 mode in all flow quantities. Higher 
order modes are most apparent in the density p and vertical 
current density j y . In the mode injected slices on the other 
hand, we observe a slight dominance of the m = 3 mode. 

To quantify the growth of non-axisymmetric modes 
within the jet, we calculate the fast Fourier transform of 
the variables on the slice in cylindrical (r, <j)) representa- 
tion. For this purpose, we re- grid the unstructured slice 
data x G [—12,12], z £ [—12,12] containing the jet spine 
to a uniform grid before transformation to the cylindrical 
coordinates r £ [0, 12], <j> £ [0, 27r]. Thereafter, the Fourier- 
transformation of the (r, 0)-plane 

N r 

/(n,m)= ^2 E f( n r,ri<t>) 



n r = -N r n (b = -N cb 



(7) 



is executed which yields the radial (n) and angular (m) 
Fourier amplitudes of the input scalar via A(n,m) = 
|/(ra,n)| 2 . To quantify fluctuations of the angular part 
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Figure 3. Slice along the jet in simulation M3D at t = 103 z = 0. 
Shown are (left to right, top to bottom): Thermal pressure (log- 
scale); co-moving density (log-scale); (toroidal) magnetic field 
strength across the plane and magnitude of four velocity u = Tv. 
Critical points are indicated in the plot of density, where we show 
the (critical) Alfen surface (blue), the light cylinder (black) and 
the fast magnetosonic surface (red) according to the definitions 
of |Camenz ind (1986). 



alone, we define the normalized cumulative Fourier ampli- 
tudes 



A{m) 



imo)i 2 n=-V n+ i 



E \f(n,m)\ 2 



(8) 



which measures the fluctuations with angular frequency m 
in relation to the squared mean of the scalar / expressed 
by |/(0, 0)| 2 . The Fourier amplitude planes of the density 
fluctuations of Figure [4] and [5] are shown in Figure [6] The 
m — 4 pollution of the unperturbed run is clearly visible, 
mode injection on the other hand can be used to get rid 
of this effect almost entirely as shown in the lower panel of 
Figure [6] 

3.2.1 Temporal Evolution 

To quantify the temporal growth of the modes, we calculate 
the Fourier amplitudes at the y — 32 slice in run L3Dm 
for various snapshots. At this altitude, the magnetic "back- 
bone" By becomes distorted at t > 80 with dominating 
m = 1 and m = 2 modes, however the amplitudes grow no 
further. The time evolution of the A(m) function is shown 
in the left panel of Figure [7] After an exponential rise where 
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Figure 6. Top: Radial (n) and angular (m) Fourier amplitudes 
of density across y = 32 at t = 100 for the unperturbed case 
L3D (left) and for the perturbed case L3Dm (right) . Bottom: 
Cumulative modes for the two cases. To guide the eye, we show 
the empirical mode-decay following the power-law m~ 3 . 
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Figure 8. Instantaneous barycenter displacements in units of 
inner disk radius and cumulative Fourier amplitudes of density p, 
all along the unperturbed jet L3D at time t = 168. The jet base 
is at left. 
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Figure 7. Left: Mode growth of B y at y = 32 in simulation 
L3Dm. After initial exponential rise, the modes tend to saturate. 
Due to grid-noise, the m = 4 mode is initially comparable to the 
dominant m = 1 mode. Right: Barycenter motion on the y = 32 
slice. 



the growth time in m — 1 is shortest, followed by the m — 4 
mode, the perturbations saturate at t ~ 80 and subsequently 
fluctuate about a mean value. Mostly, the amplitudes are or- 
dered according to A(m) > A(m +1) although the m — 2 
occasionally surpasses the m—\ contribution. 



3.2.2 Spatial Evolution 

The clearest indicator of the kink instability can be observed 
in the deflection of the jet barycenter. For this purpose we 
define the barycenter f = 



\J x 1 + z 1 of the quantity Q 

dx dz _ _ f z Q dx dz 



f Q dx dz f Q dx dz ^ 

in analogy with Mignone et al. ( 2010| ). For the density- 
displacement r cm , we define Q cm = X^p, where the x is 
computed from the tracer scalars to pick out the jet contri- 
bution alone. This quantity is also shown against simulation 



time in the right panel of figure [7] For the current displace- 
ment Tj y> o we set Qj y >o — xjy taking only the positive 
values of the current density jy . Finally, the motion of the 
magnetic flux is defined via Qb v >o = X^>y , taking into ac- 
count only the positive flux . Instantaneous barycenter 
position and mode population along the jet is shown in Fig- 
ure [8] for run L3D and in Figure [5] for the mode-injected run 
L3Dm. Let's first focus on run L3D. The barycenter dis- 
placement is small compared to the inner disk radius and 
even tends to decrease along the jet. Only near the jet-head, 
at about y ~ 80 a significant displacement can be observed. 
The modes are dominated by ubiquitous m = 4 noise, which 
seems to suppress all other fluctuations starting at height 
y — 20. The m — 4 dominance prevails all the way to the 
jet head. 

The behavior of run L3D is different, here, the kink 
mode surpasses the m = 4 at y ~ 10. The angular fluctu- 
ations saturate around y — 20. Until y = 60, the displace- 
ments in current and magnetic flux stay roughly constant. 
This hints to a self-stabilization of the jet formation region. 
The kink mode starts to rise again towards the jet head, 
accompanied with a notable barycenter displacement. We 
note that the magnetic field configuration near the jet head 
is strongly toroidally dominated as the monopolar initial 
field configuration is overtaken by the jet (see e.g. Figure 
[2|. This toroidal dominance could yield an explanation for 
the stronger growth of the kink mode at the jet head. To 
quantify this further, we introduce the co-moving magnetic 
pitch defined as 



(10) 



which plays a major role in the stability of current carrying 
plasmas in the laboratory (e .g. |Bateman 1978|) and astro- 
physical jets (e.g. jAppl et al.||2000||Lery et al.||2000| ). Small 
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Figure 9. As Figure [8] but for run L3Dm at time t = 160. The 
mode injection at the base (y=0) avoids artificial m = 4 domi- 



values of P/r < 1 and thus toroidally dominated configu- 
rations are particularly susceptible to the kink instability. 
The co-moving fields are obtained by applying the projec- 
tion B' a = upF* afi , where F* a(3 denotes the dual Faraday 
tensor and up the co- variant Four- velocity as customary. 
The radius r and the fields , B p involved in the definition 
of the pitch are well-defined however only in axi-symmetry. 
For small perturbations from the cylindrical shape, we can 
re-orient the symmetry axis to the magnetic backbone at 
the position (x, z) and define effective values for f, and 
Br with respect to the new origin. To define the magnetic 
backbone position using equations (|9j, we apply the kernel 
Qbb — xBy which reliably finds the peak of magnetic flux. 
To investigate also effects due to the electric field, we de- 
fine an effective light surface via the comparison of the field 
strengths 

(id 



Bp 

where x = 1 marks the light surface and Bp — yj B^ + B% ~ 
By is the poloidal field strength with respect to the mag- 
netic backbone. Since the jet is well collimated B y ^> B f , 
the location of the backbone is in fact only of secondary 
importance for the definition of the light surface and we ob- 
tain similar results when considering only the vertical field 



By in Equation (11). This light surface location is indicated 
with the gray contour in the top panel of figure [To] where 
the pitch profiles across the jet is shown. To visualize the 
Lorentz force across the flow, we show force vectors of the 
electric field p e E = (V • E) E and of the classical Lorentz 
force j x B where we neglected the displacement- current for 
simplicity. 

The axial flow is outside of the light surface, then fol- 
lows a region of super-luminal field line rotation and a third 
region outside of the main jet where field lines rotate again 
sub- luminal. The magnetic backbone is markedly seen in 



-10 - 




Figure 10. Top: Co-moving pitch and effective light surface with 
respect to the center of magnetic flux (shown as grey contours) 
through the surfaces y = 60 for run L3Dm at t = 160. Black 
arrows indicate the direction of the electric force p e E and white 
arrows show the direction of the Lorentz force (V x B) x B. 
Both force vectors are projected onto the image-plane. The axis is 
marked by a black "+" and the center of magnetic flux is shown as 
gray "+" . Bottom: Corresponding magnetic flux B y with velocity 
(white) and magnetic field vectors (black) in the plane. 



the pitch and we note that electric forces exhibit a de- 
collimating contribution at the central core while they tend 
to collimate near the outer light cylinder. In most regions, 
the additional electric component is opposed to the classi- 
cal Lorentz force. As a general trend, we obtain a radially 
decreasing pitch profile from the backbone to its minimal 
value P/r < 1, from where the pitch increases again to the 
boundary of the jet. Due to the m=l motion, a spiral pat- 
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Figure 11. Field lines of the unperturbed (top) and perturbed 
(bottom) simulations at times t = 103 and t = 194, respectively. 
To guide the eye along the bent jet, pressure isocontours are added 
in the figures. The perturbed run shows a wider magnetic back- 
bone and decreased toroidal field. 



tern with increased pitch is shed of the magnetic backbone. 
This 'smearing' leads to a diminishing of the toroidally dom- 
inated regions. These regions of low pitch coincide with the 
locations of super-luminal field-line rotation and the stabi- 
lizing effect of electric forces becomes apparent as they tend 
to counter-act the magnetic contribution. 



3.3 Jet-Cloud interaction 

We now study how the jet reacts to external perturba- 
tions exerted by a clumpy ambient medium. As the jet head 
punches its way through the in- homogeneous environment, 
also the upstream flow becomes deflected. Figure [TT] shows 
a rendering of the simulation M3Dmd with mode injection 
and clumpy medium in comparison to the unperturbed case 
M3D. Field lines are colored according to the magnetic pitch. 
The clouds with a density contrast of 100:1 represent a 
strong perturbation and we find the jet heavily distorted. 
Accordingly, the motion of the jet bary center is increased 
and we find a strong dominance of the m — 1 mode along 
the whole jet (Figure 12). Due to the increased density of 




the external medium, the jet propagation speed is reduced 
and the bow-shock is much wider (compare with Figure 
[2|. We therefore compare the jet morphology between runs 



0.001 



Figure 12. Barycenter motion and low order angular modes in 
density of simulation M3Dmd at t = 194. 



M3D and M3Dmd at times of roughly equal jet propagation 
length. The precession of the magnetic backbone against the 
toroidal magnetic field direction tends to "smear" out the 
high-pitched axial region and the tightly wound helix is ef- 
fectively unwound. This represents an efficient mechanism 
of jet self-stabilization also noted by |Ouyed et al.| ( f2QQ3| ) . In 
their study the effect was described as follows: "The appear- 
ance of the \m\ = 1 modes pumps energy into the poloidal 
magnetic field, causing the jet Alfven Mach number to fall 
below unity and stabilize the jet" . From our simulations we 
come to a similar conclusion, in addition, we note that the 
"unwinding" of the helical field also decreases the field-line 
rotation Q, and thus also the influence of electric fields. We 
show the increase in magnetic pitch compared to the unper- 
turbed simulation L3D in Figure [13] Regions of low pitch 
are reduced to a thin sheet where the effective light cylinder 
is found. In this regime, electric stabilization can not be of 
importance since the typical value of the light-cylinder x is 
only of order unity or below. 

It is interesting to visualize the return current in the 
heavily perturbed run M3Dmd. The structure of the current 
density across the jet is shown in figure |14| Fine layered 
current sheets form on small scales that could in principle 
dissipate and give rise to particle acceleration within the jet. 
It is only due to the high resolution that this effect can be 
observed. 



4 CONCLUSIONS 

We have presented first results of high-resolution 3D sim- 
ulations of relativist ic jet formation from magnetospheres 
in Keplerian rotation. When the flow is perturbed by non- 
axisymmetric internal perturbations of the accretion disk 
corona, the modes first grow exponentially at the base, ap- 
proach saturation along the jet and grow again towards the 
jet head. The m — 1 kink is the dominant mode of departure 
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Figure 13. As in Figure |10| comparison of the pitch at y = 
32 in runs M3Dmd (left) and L3Dm (right) at times t = 194 
respectively t = 160. The motion of the backbone in the heavily 
perturbed simulation increases the poloidal field on account of 
the unwinding of the helical field. 
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Figure 14. Out of plane current density at y = 32, t = 194 for the 
run with external perturbation M3Dmd. A layered filamentary 
structure of alternating current directions develops. 



from axisymmetry. At a given height above the accretion 
disk, the temporal evolution of the modes was considered. 
Also here we find a saturation of perturbations before a no- 
table dissipation or even disruption is encountered. As an 
aside, we also performed simulations where the only mea- 
sure of perturbation is the ubiquitous discretization noise 
of the grid. In result, the modes are dominated by multi- 
ples of m — 4 which grow along the flow at the expense of 
other modes and we observe virtually no motion of the jet 
barycenter. 

To further investigate the stability of the jet structure, 
we considered the co-moving magnetic pitch. As in the ax- 
isymmetric case, a stabilizing backbone of high pitch P ^> 1 
develops, surrounded by an intermediate, toroidally dom- 
inated region P < 1 and an outer high-pitch region at 
the border of the jet. The locations of super-luminal field 
line rotation (the light surface) approximately coincide with 
the low-pitch region. Forces due to the electric field p e E 
oppose the classical magnetic Lorentz force (V x B) x B 



which could thus add to jet stabilization in the relativistic 
case. 

In order to study external perturbations, we initial- 
ized the domain with a static clumpy ambient medium fol- 
lowing a power-law spectrum in Fourier space. While this 
should not be mistaken for a fully developed MHD-turbulent 
medium, it allows us to investigate perturbations "external" 
to the jet as a general scenario. The maximum amplitude in 
cloud density to the jet density was chosen as 100 : 1 which 
thus represent a strong perturbation to the emerging jet. 

In result, the jet funnels its way through the path 
of least resistance which leads to large departures of the 
barycenter from the axis and dominating m = 1 modes. 
When compared to the unperturbed case, the magnetic pitch 
is largely increased which can be interpreted as a mechanism 
of jet self-stabilization. Due to the precession of the mag- 
netic backbone against the toroidal magnetic field direction, 
the helical structure tends to be "unwound" leading to an 
increase of poloidal field which also reduces the amount of 
field-line rotation. The external perturbation and accompa- 
nying motion of the magnetic backbone also gives rise to 
filamentary small-scale structure visible in the return cur- 
rents. The accompanying dissipation could facilitate particle 
acceleration within the jet. Further investigation with high- 
resolution 3D simulations is needed to quantify this effect. 
Another point of interest is the validity of the axisymmetric 
assumption in the formation of jets in general. We find that 
due to the early saturation of non-axisymmetric instabili- 
ties, 2D results on the inner scales of jet acceleration and 
collimation (e.g. PF I, PF II) appear largely unaffected by 
non-axisymmetric effects. Detailed comparisons of the dy- 
namics in two- and three- dimensions will be provided in a 
forthcoming paper. 
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